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Abstract. The phenomena of stochasticity in biochemical processes have been intrigu- 
ing life scientists for the past few decades. We now know that living cells take advantage 
of stochasticity in some cases and counteract stochastic effects in others. The source of 
intrinsic stochasticity in biomolecular systems are random timings of individual reactions, 
which cumulatively drive the variability in outputs of such systems. Despite the acknowl- 
edged relevance of stochasticity in the functioning of living cells no rigorous method have 
been proposed to precisely identify sources of variability. In this paper we propose a 
novel methodology that allows us to calculate contributions of individual reactions into 
the variability of a system's output. We demonstrate that some reactions have dramati- 
cally different effects on noise than others. Surprisingly, in the class of open conversion 
systems that serve as an approximate model of signal transduction, the degradation of 
an output contributes half of the total noise. We also demonstrate the importance of 
degradation in other relevant systems and propose a degradation feedback control mecha- 
nism that has the capability of an effective noise suppression. Application of our method 
to some well studied biochemical systems such as: gene expression, Michaelis-Menten 
enzyme kinetics, and the p53 system indicates that our methodology reveals an unprece- 
dented insight into the origins of variability in biochemical systems. For many systems 
an analytical decomposition is not available; therefore the method has been implemented 
as a Matlab package and is available from the authors upon request. 



Living cells need to constantly adapt to their changing environment. They achieve 
this through finely honed decision making and stress response machineries that regulate 
and orchestrate the physiological adaptation to new conditions. In all studied genomes 
a large number of proteins have as their primary function the transfer and processing 
of such information. Such proteins are linked through a host of different mechanisms 
into biochemical circuits that perform a variety of information processing tasks including 
storage, amplification, integration of and marshalling the response to environmental and 
physiological signals [I] . The functioning of these information processing networks depends 
on thermal or probabilistic encounters between molecules, resulting in a distortion of a 
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transferred information that is best understood as noise. Each reaction in the information 
processing machinery leads to an inevitable loss of information [2]. Therefore, cell functions 
do not only rely on the necessity to make "good" decisions, but also on appropriate ways to 
cope with the uncertainties arising from the "noisy" signal transmission. To deal with the 
latter type of difficulty, we believe, evolution equipped cells with reliable signal transduction 
systems by using less noisy reactions or reaction configurations, where needed [3]. The 
question, however, which reactions, molecular species or parts of a network contribute 
most of the variability of a system or are responsible for most of the information loss has 
not gained much attention, except for some models of gene expression [U EJ EE [7] . Origins 
of stochasticity in biochemical systems have been discussed and argued over at length, but 
a unified and robust mathematical framework to study this problem has been lacking. 

Here we present a novel and general method to calculate contributions of individual 
reactions to the total variability of outputs of biochemical systems. This enables us to 
identify the origins of cell-to-cell variability in dynamical biochemical systems. We derive 
a modified fluctuation-dissipation theorem which enables us to determine how much of the 
total variance results from each of the system's reactions. We then obtain some unexpected 
but general rules governing the noise characteristics of biochemical systems. In particular, 
we shall show that in an arbitrary system with a sufficiently large number of molecules, 
degradation of the output (e.g. a reporter protein or a transcription factor) contributes to 
the total variance of the system half the of the output's mean; for the important class of 
open conversion systems exactly half of the variance derives from the degradation step of 
the output signal. These results demonstrate that some reactions may be responsible for 
higher information loss than others; but our results also reveal that cells have the option 
of optimising biochemical network structures in order to avoid the most noisy reactions if 
necessary. Based on these results we propose a mechanism of controlled protein degradation 
based on a positive feedback that allows an arbitrary noise reduction resulting from the 
protein degradation. 

Below we first introduce the general framework for modelling chemical reactions and 
derive a new method to decompose the noise in a biochemical system into contributions 
from different individual reactions. Furthermore, two general properties governing noise 
are presented. Finally, we use biological examples of signal transduction systems to provide 
novel insights into the origins of variability. In particular, we decompose the variance of the 
outputs of linear transduction cascades and Michaelis-Menten enzyme kinetics. In addition, 
for the oscillatory p53 system we show that stochasticity in p53 protein degradation is 
responsible for the variability in the oscillations' periodicity. 

1. Noise Decomposition 

We consider a general system of n chemical species [8] inside a fixed volume with 
x = (xi, . . . , x n ) T denoting the number of molecules. The stoichiometric matrix S = 
{sij}i=i.2...n; j=i,2...r describes changes in the population size due to r different reaction 
channels, where each Sij describes the change in the number of molecules of type i from Xi 
to Xi + caused by an event of type j. The probability that event j occurs during the 

2 



time interval [t,t + dt) equals fj(x,t)dt, where the fj(x,t) are called the transition rates. 
This specification leads to a Poisson birth and death process described by the following 
stochastic equations |9j 



(1) x(t) = x(0) + j2s. j N j (f f j (x,s)ds), 

where the Nj are independent Poisson processes with rate 1. 

In order to define the contribution of the jth reaction (j = l,...,r) to the variability 
of x(t) we first define as the expectation of x(t) conditioned on the processes 

N%(t), iVj_i(i), Nj+i(t), N r (t), so that {x(t))fj\ is a random variable where timings 
of reaction j have been averaged over all possible times, keeping all other reactions fixed. 
Therefore x(t) — (x(t))ij is a random variable representing the difference between the 
native process x(t) and a time-averaged process of the jth. reaction. Now the contribution 
of the jth reaction to the total variability of x(t) is 

(2) E«)(i) ee ((x(t) - (x(t)) H )(x(t) - (x(t)) H f), 

where (...) denotes the temporal average over all r reactions. This definition is similar 
to the one proposed in [7J to quantify the contributions of promoter states and mRNA 
fluctuations, respectively, to the protein level variability. 

In general, it is difficult to calculate or study properties of equation Q using a Poisson 
birth and death process framework . Here instead we use the linear noise approximation 
(LNA), which allows us to model stochastic systems using Wiener processes driven (Wj(s)) 
instead of Poisson processes driven (Nj(s)) stochastic differential equations [El [TOJ HH 
I12j . The LNA is valid if the number of interacting molecules is sufficiently large |12j and 
decomposes the system's state, x(t), into a deterministic part ip(t) and a stochastic part 

m 

(3) x(t) =<p(t) + £(t); 

here ip{t) and £(t) are described by the deterministic and stochastic differential equations 

(4) <p(t) = p(0)+ [ SF(<p(s),s)ds 



(5) + J2 s -i w 3 y/fM*)>s)*>) , 

respectively, and their coefficients are given by the following formulae 

(6) F( l p,t) = (f 1 (<p,t),...,f r ( V ,t)) T 

(7) { A M}ik = J2 s ^ k 
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The LNA presents a simple way to compute contributions and here we demonstrate 

how the total variance can be decomposed into the sum of individual contributions. We 
first write the explicit solution for the process £ as 

(8) m = $(0, i)£(0) + / $( s , t)E(<p, s)dW(s), 

Jo 



and E{ip,t) = S diag y y/ fi(<p, t)), \J (f r (<P, t)j , where t) is the fundamental matrix 
of the non-autonomous system of ordinary differential equation^] 

fill's t) 

(9) a L±±A =A (t)<t>(s,t), *(8,8) = I. 
Now it is straightforward to verify that 

(10) <£(*)>[_,•= *(0,t)£(0)+ f $(s,t)E^(s)dW(s), 



where (t) = S™E(t), E®(t) = S®E(t), = {i~ i) }/=i,..,n,fc=i,...,r, 

{s}fc }i=l,..,n,fc=l,...,r and 

(-i) _ / for A; = j m _ / si k for A; = j 



Slk l for A; 7^ j S/fc 1 for A; / j 



From (JSj) and (10) we have 

(11) m- < m > hJ = f t Hs,t)E^( S )dw(s). 

Jo 

With x(t) - {x(t))\_j = - (£(i))|_j and the time derivative of ((£(*) - (£(i))|_j)(£0) - 
(£M)|-j) T > we obtain for S^(i), 

TV"! (?) 

(12) — — = A(i)E&') + £ (j 'U(t) T + £> (i) (t), 
with 

(13) DW(t)=EV\t)(E®(t)) T . 

This is, of course, analogous to the fluctuation dissipation theorem, with the exception that 
the diffusion matrix contains zeros at all entries not corresponding to the jth reaction. Now 
the fact that the total variance 

E(t) = <(*(*) - <s(t)»(x(t) - (x(t))f) 

can be represented as the sum of individual contributions 

(14) £(i) = £«(*)+ ... +£ (r) (t), 



^In order to simplify notation, form now on we will write A(t) and E(i) instead of A(ip,t) and E(ip,t), 
respectively. 
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results directly from the decomposition of the diffusion matrix D(t) = $^ = i ■D^(i) and 
the linearity of the equation for S(t), given by the standard fluctuation dissipation theorem 

(15) ^ = A(t)Z + ZA(t) T + D(t). 

2. Mathematical Results 



With the decomposition (14) it is in principle possible to detect reactions that make 
large contributions to the output variability of biochemical reactions. But even simple 
systems, for which analytic formulae exist, usually have complicated noise structures; we 
can nevertheless prove two general propositions that assign surprisingly substantial contri- 
butions to the degradation of an output signal. We formulate them in this section (proofs 
are in the Appendix) and illustrate them further below. 

2.1. Proposition 1. Consider a general system such as described at the beginning of the 
Noise Decomposition section. In addition, assume that the deterministic component of 
the system ip has a unique and stable stationary state (all eigenvalues of matrix A have 
negative real parts). If x n is an output of this system being produced at rate Yl,i^n kf x i an d 
degraded in the rth reaction at rate k~x n then the contribution of the output's degradation 
is equal to the half of its mean; more specifically, 

(16) [EW] nn = ±(x n ). 

2.2. Proposition 2 (open conversion system). Now consider again a general system 
but assume that reaction rates are of mass action form and that only three types of reactions 
are allowed: 

k + 

production from source:: — x% 

k~ 

degradation:: X{ — 

conversion:: X{ %i '> Xj. 

To satisfy the openness assumption each species can be created and degraded either directly 
or indirectly (via a series of preceding conversion reactions) . As in Proposition 1 let x n be 
the output of the system. Under this assumption the degradation of the output contributes 
exactly half of the total variance of the system's output, 

(17) [s«l =hn nn , 

L J nn Z 

where r is again the index of the output's degradation reaction. 

2.3. Remarks on propositions. Proposition 2 can be understood as a balance between 
production and degradation reactions. If we consider all reactions except rth as producing 
species x n , then production and degradation contribute the same amount of noise. Usually, 
however, there is more than one production reaction, and therefore it is more convenient 
to interpret this result as the contribution of a single reaction. 
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Both propositions indicate that a substantial part of noise is a result of the signal's degra- 
dation reaction. This observation is particularly interesting as open conversion systems are 
often good approximations to important biochemical reactions such as Michaelis-Menten 
enzymatic conversions or linear signal transduction pathways, which we discuss below; it 
suggests also that controlled degradation is an effective mechanism to decrease the overall 
contribution to variation in protein levels. 

3. Biological Examples 

In this section we demonstrate how our methods can be used to study origins of noise 
in biochemical systems and present consequences of the above propositions. 

3.1. Simple birth and death process. For illustrative purposes we start with a basic 
birth and death process, where molecules x arrive at rate k + and degrade at rate k~x. In 
the LNA this process can be expressed by the following stochastic differential equation 

(18) dx = (k + - k~x{t))dt + \/fc+dWi + ^k-{x{t))dW 2 . 

V v ' v ' 

birth noise death noise 

The stationary distribution of this system is Poisson with the mean (x) = fp. Therefore in 
the stationary state death events occur at rate k~ (x) which must be equal the birth rate k + . 
The noise terms in the above equation are equal at stationarity indicating that contributions 



of birth and death reactions are equal. Using formula (12) it is straightforward to verify 
that the decomposition 

(19) £= \{x) + \{x) 

birth noise death noise 

holds indeed. 

3.2. Linear signal transduction cascades. Here we first focus on conversion path- 
ways, where molecules x\ are born at a certain rate kf and information is transmitted to 
molecules x n using (possibly reversible) conversion reactions. All molecules can degrade 
with arbitrary first order rates, and we assume that the final product x n is not the sub- 
strate of any conversion reaction. Under these assumptions, regardless of parameter values, 
length of the pathway, and degradation of intermediates, the degradation of the output x n , 
rather unexpectedly, contributes exactly half of the variance of the signalling output, x n . 
On the other hand, if information is transmitted not using conversion reactions, but each 

species catalyses the creation of the subsequent species, i.e. x\ - 1+1 '> Xi + Xi+i, then 
according to the Proposition 1, the contribution of the degradation of x n to the output's 
variance equals half of its mean regardless of the cascade length and parameter values. 

In order to study contributions of intermediate reactions, we considered cascades of the 
length three (see Figure [TJ and decomposed variances numerically for both conversion and 
catalytic cases. Decompositions are presented in Figure [2] For conversion cascades we 
used two parameter sets corresponding to fast and slow conversions to demonstrate that 
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increasing the rates of intermediate steps decreases their contribution as might intuitively 
be predicted (top panel of Figure [2]). For catalytic cascades increasing the rates of reactions 
at step two and three (going from slow to fast) increases the contribution of reactions at the 
bottom of the cascades. For slow dynamics during steps two and three, the relatively fast 
fluctuations at the start are filtered out (low pass filtering) . If the dynamics of all species 
occur at similar times scales then these fluctuation can efficiently propagate downstream 
(bottom panel of Figure [2]). 

3.3. Michaelis-Menten kinetics. Michaelis-Menten enzyme kinetics are fundamental 
examples of biochemical reactions; substrate molecules (S) bind reversibly to enzyme 
molecules (E) with the forward rate constant ko and the backward rate constant k\ to 
form a complex (C), which then falls apart into the enzyme and a product (P) at rate k 2 - 

k S-E k2 

S + E C — > E + P. To ensure existence of the steady state, we assume that sub- 

fci 

strate molecules arrive at rate kb and are degraded at rate kd- At the unique steady state, 
if the concentration of an enzyme is large compared to the substrate, the system is well 

k S-E k2 

approximated by a set of mono- molecular reactions: S C — > P. Our theory therefore 

fci 

predicts that half of the noise in Michaelis-Menten kinetics operating in its linear range 
(abundant enzyme) is generated by degradation of the product molecules. In Figure [3] we 
have calculated the contributions of each of the four reactions to the variance of all four 
species for the full model (without themonomolecular approximation). The contribution 
of the output degradation is as predicted by the monomolecular approximation. 

3.4. Gene expression. The canonical example of a linear catalytic pathway is the gene 
expression process, which can simply be viewed as the production of RNA (27) from source 
at the rate kf, and production of protein (x 2 ) in a catalytic reaction at rate k 2 x\ together 
with first order degradation of both species at rates k±x\, k^ x 2 . Proposition 1 states that 
the part of the variance resulting from the protein degradation equals half of the protein 



mean. Moreover, formula (12) allows us to derive the complete decomposition 




(20) V(x 2 ) = -(x 2 ) + -(x 2 ) + + lk ^) 




2 7 r + 7 p ' 2 7 r + 7 p 

prot. degradation translation mRNA degradation transcription 

A similar decomposition was first presented by [61 [5], however, only into contributions 
resulting from fluctuating species and without theoretical explanation; the latter was pro- 
vided later by [7j. Below we investigate two extensions of the above model. First, we 
assume that the promoter can fluctuate between "on" and "off" states (similarly to P2E]) 
and calculate contributions for different time-scales of these fluctuations. Second, we as- 
sume that the protein is a fluorophore that undergoes a two-step maturation process before 
it becomes visible (folding and joint cyclization with oxidation). Figure [4] presents contri- 
butions for fast and slow promoter kinetics, showing that fast fluctuations are effectively 
filtered out (contributing 10 %) but remain a substantial contributor when they are slow 
(contributing 40 %). 
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Variability in gene expression is often measured by means of fluorescent proteins that 
undergo maturation before becoming visible for detection techniques; but the process of 
maturation |13| itself is subject to stochastic effects, and can thus contribute significantly 
to the observed variability. We used typical parameters |14[ [13] for fast and slow maturing 
fluorescent proteins and found that maturation contributed 4 and 25%, respectively (Figure 
[5]) to the overall variability; here our method allows for the rigorous quantification compared 
to the previous qualitative observations of |15|. [16] . 

3.4.1. Out- of- steady state systems: p53 model. In order to demonstrate that our method- 
ology and predictions are valid also for more general out-of-steady state models we next 
focus on the p53 regulatory system which incorporates a feedback loop between the tumour 
suppressor p53 and the oncogene Mdm2, and is involved in the regulation of cell cycle and 
the response to DNA damage. We consider the model introduced in [T7] and later analysed 
in [18] that reduces the system to three molecular species, p53, mdm2 precursor and mdm2, 
denoted here by <p p , (p yo and <p y . The model incorporates six reactions (i = 1, . . . , 6), and 
its deterministic version can be written in the form of the following ordinary differential 
equations 

(21) (p p = h - k 2 (f p - facfy — ^— - 

(fp + H 

(22) (p yo = k A ip p - k 5 (p yo 

(23) (p y = k5(p yo - kefy. 



Using equation ( 12 ) we calculated the contributions of all of the six reactions present in the 
model. Results are presented in Figure [6} The contributions of each reaction oscillates over 
an approximately constant range over the course considered here, except for the reaction 
of p53 degradation: the contribution of the p53 degradation accumulates over time. This 
observation is consistent with our propositions which predict for the steady-state models 
that the contribution of output signal degradation is significant, and generally overrides 
the contributions of other reactions. 

In [19] the authors studied the applicability of the LNA for the analysis of oscillatory 
models and found that the approximation fails in this type of models when long time 
periods are considered. In these cases the total variance diverges to infinity with time. 
They also decomposed the total variance into an oscillatory and a diverging part; the 
diverging part corresponds to the variability in the oscillation's period. The diverging part 
in Figure [6] results from the p53 degradation, indicating that stochasticity of this reaction 
is responsible for the variability of the period in this system. 

4. Reducing the contribution of degradation noise 

The propositions exemplified by the above models demonstrate that the noise resulting 
from degradation of an output is a substantial source of variability in biochemical systems. 
In a general system it contributes a fraction of the variance equal to the half of the mean 
output. Here we show how controlling the degradation can reduce this contribution. We 
consider a simple birth and death process, but the same mechanism is valid for general 



system considered in the Proposition 1. If births and deaths occur at the state dependent 
rates f(x) and g(x), respectively, such a system is described by 

(24) dx = (f(x) - g(x))dt + y / f((x))dW 1 + y/g{(x))dW 2 • 



birth noise death noise 



At equilibrium, we have f({x)) = g((x)), and using formula (12) we can calculate that 
contribution of degradation is 

/r>r\ y(death) _ \9{{ x )) 

K ' df((x)) dg((x)) - 

dx dx 

For f(x) = k, g(x) = ^yx this reduces, of course, to the previously discussed case 

S (death) = 

\{x). Nevertheless, if functions / and g are of Hill and Michaelis-Menten type, respec- 
tively, then the contribution is no longer directly related to the mean and can be reduced 
to an arbitrary low level according to the above formula. Any decrease in the contribution 
compared to the mean is achieved by the reduction of the variance resulting from the re- 
duced flux through both reactions and autoregulatory control effects. Recent experimental 
work on protein degradation in vitro provides the evidence that degradation indeed ex- 
hibits Michaelis-Menten type kinetics [20], instead of the linear first order kinetics that are 
usually used to model degradation. The effect of autoregulation is depicted in Figure [7| 



5. Discussion 

The noise decomposition method introduced here allows us to investigate in detail where 
and how noise enters biochemical processes, is propagated through reaction systems, and 
affects cellular decision making processes. We have shown analytically that in a wide class 
of systems the signal degradation reaction contributes half of the noise in a system's output 
regardless of parameter values. We have also carried out numerical study for a system that 
never reaches steady-state and confirmed the surprisingly important role of degradation in 
controlling stochasticity of the p53 system, as well as in a range of generic systems biology 
models. Quite generally, the ability to dissect noise propagation through biological systems 
does enable researchers better to understand the role of noise in function (and evolution), 
and will also enable synthetic biologists to either harness or dampen the effects of noise in 
molecular signalling and response networks. 

One of the central results that we found and report here is the crucial role of degradation 
of the signal on the overall noise levels. The relevance of degradation has not been studied 
before in the context of stochastic biochemical dynamics and to our knowledge our study is 
the first that draws attention to the importance of degradation. This is particularly impor- 
tant as this may indicate new therapeutic targets: in humans and other other sequenced 
organisms, certainly, the repertoire of proteins involved in protein degradation, in partic- 
ular ubiquitin-ligases, is as rich and diverse as the repertoire of proteins regulating their 
activation, the kinases (and phosphatases). Thus targeting the degradation of proteins 
appears as important to biological systems as protein activation and offers an attractive 
and broad range of new potential therapeutic targets. 
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Appendix A. Proof of Proposition 1 

1) Interactions of x n with other species imply that a%j = for i = 1, ...,n — 1 and j = n. 
Thus matrix A can be written as 

A _ ( A 



2) Formula (13) implies that all elements of matrix are equal to except D 
kh \ x n), therefore has the form 

D« = ( ° ° 

3) It is straightforward to verify that the n x n matrix 

< 26 > *" - ( o s<:„> 

satisfies the equation 

(27) + ^A T + DM = o 



(r) 
mi 



Appendix B. Proof of Proposition 2 

In an open conversion system S nn = (x n ) (see e.g. [21]), therefore from Proposition 1 we 
have Yj nn = 2Yj nn . 
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Figure 1. Three step conversion (A) and catalytic (B) pathways. Rates fj 
for reactions j = 1, 6 have the form: f\ = kf , /2 = k^xi, /a = k^X2, f& = 

^4 x Xi = ^5 %5i ffi = kg Xq. 



conversion - slow conversion - fast 




Figure 2. Variance contributions of the different reactions in the three 
steps linear conversion (top) and catalytic cascades (bottom). Parameters 
for each case are given below. Conversion slow: kf = 50, /cj = l,k% = 
1,^4 = 1,^5 = I, kg = 1; conversion fast: kf = 50, k^ = 10, k% = 
10, fcj = 1, fcg = = 1; catalytic slow: A;^ = 50, k J = 0.1, k% = 

0.1, k^ = 1, fcj" = 0.1, k$ = 0.1; catalytic fast: fc^ = 50, = 10, fog = 
10, fcj = 1, fcrj = 10, fc^ = 10. All rates are per hour. 
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substrate enzyme 




Figure 3. Variance decomposition for the four species of the Michaelis- 
Menten kinetics model. The following parameters were used: ko = 0.1, k\ = 
50, k% = 1, kf, = 20, kd = 0.1. All rates are per hour. 




29% 
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Figure 4. Variance contribution of the reactions involved in gene expres- 
sion with fluctuating promoter states with transcription rate modelled as 
k~l = i+yjlj , arid y being the stationary solution of dy = (b — ^y)dt + 
yj2rfdW . The following parameters were used: V = 100, H = 50, k^ = 2 
k^ = 0.44, &2 = 0.55. For fast fluctuations we used b = 50,7 = 1) an d for 
slow b = 0.5, 7 = 0.01. All rates are per hour. 
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Figure 5. Variance contributions of the reactions involved in fluorescent 
protein maturation model for slow (left) and fast (right) maturing proteins. 
The following parameters were used: kf = 50, k± = 0.44, kf = 2, k^ = 
0.55. For slow maturation (average maturation time approx. 5 h.) we 
assumed folding and maturation rates to be kf = 0.2 and k m = 1.3, respec- 
tively. For fast maturation (average maturation time approx. 0.5 h) we set 
kf = 2.48, k m = 13.6. All rates are per hour. 




Figure 6. Trajectories (top left) and decomposition of variance into the 
contributions corresponding to each of the 6 reactions for mdm2 (top right), 
mdm2 RNA (bottom left) and p53 (bottom right). We used the index i of 
the parameters ki in equations ( 
We used model and parameters published in [17]. Figure demonstrates 
accumulation of the noise contributed by degradation of p53 (6th reaction) . 
The following parameters were used: k\ = 90, ki = 0.002, kj, = 1.7, H = 
0.01, ki = 1.1, k§ = 0.8, k% = 0.8. All rates are per hour. 



21p3) to denote reaction they describe. 
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Figure 7. The effect of noise reduction resulting from regulated degrada- 
tion. Each colour describes production rate f(x) (dotted line) , degradation 
rate (dashed line) and density of stationary distribution of x (solid line). 
For blue: f(x) = 1/2, g(x) = 0.025x; for black: f(x) = 1/(1 + (x/20) 5 ), 
g{x) = (x/20) 5 /(l + (z/20) 5 ); for green: f(x) = 1/(1 + (x/19) 50 ), g(x) = 
(x/21) 50 / (l + (x/22) 50 ). Plotted densities are kernel density estimates based 
on 10000 independent stationary samples generated using Gillespie's algo- 
rithm [22]. 
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